import sys
import MySQLdb
import math
import time

conn = MySQLdb.connect (host = "saskatoon",
                           user = "dudley_omics",
                           passwd = "dud1ey",
                           db = "omics")
cursor = conn.cursor ()
merged_snps_hash = {}
type_hash = {}
type_hash["synonymous"] = "Synonymous"
type_hash["nonsynonymous"] = "Missense"
type_hash["aa_deletion"] = "Indel"
type_hash["frameshift"] = "Nonsense"
type_hash["INTERGENIC"] = "Unknown"
type_hash["aa_insertion"] = "Indel"
mutation_repo = ''

def doMerge(outpath, ycr, changesfile1, changesfile2, inputhash):
	global mutation_repo
	mutation_repo = outpath
	print "calling merge snps %s %s save to %s" %(changesfile1, changesfile2, outpath)
	return merge_snps(ycr, changesfile1, changesfile2,inputhash)

def getMutationType(cgskey, chr, refstart, refend):
        cursor.execute ("SELECT ctype FROM omics.cgs_changes_classifications WHERE cgs_key = '" + cgskey + "' and chrom = '" + chr + "' and ref_start = " + refstart + " and ref_end = " + refend )
        #'SELECT VERSION()")
        row = cursor.fetchone ()
	if (row == None):
		print "Not in DB: %s %s %s" %(cgskey, chr, refstart) 
		return "unknown"
	else:
        	return row[0]

def merge_snps(crossing, ra_changes_file1, ra_changes_file2,outhash):
	global mutation_repo
	print "Begin merge %s" % time.strftime("%c")
	#/solexa/dudley/genome_assembly/YO490/YSEQ6/unprocessed_reads/alignment_1/YCGS5/RA_changes.txt
	pathspl = ra_changes_file1.split("/")
	pathspl2 = ra_changes_file2.split("/")
	cgskey1 = pathspl[len(pathspl)-2]
	merged_file = open(mutation_repo + "/" + crossing + "_" + pathspl[8] + "_" + pathspl2[8] +  "_merged.mut", "w")
	merged_file.write("chr\tstart\tend\tsample\ttype\tnotes\n")
	print crossing + " changes1 " + ra_changes_file1 + " changes2 " + ra_changes_file2
	#loop thru file one and put first five columns as key, parse strain name for sample type
	file1 = open(ra_changes_file1, "r")
	strain = pathspl[4] + ":" + pathspl[8]
	outhash["YCGSKey1"] = cgskey1 + "_" + pathspl[5] + "_" + pathspl[4]
	for line in file1:
		tokens = line.split("\t")
		mykey = ":".join(tokens[0:5])
		merged_snps_hash[mykey] = strain
        cgskey2 = pathspl2[len(pathspl)-2]
	strain = pathspl2[4] + ":" + pathspl2[8]
	outhash["YCGSKeys"] = cgskey1 + "_" + cgskey2
	outhash["YCGSKey2"] = cgskey2 + "_" + pathspl2[5] + "_" + pathspl2[4]
	file1.close()
	file2 = open(ra_changes_file2, "r")
	dupe = 0
	for line in file2:
		tokens = line.split("\t")
                mykey = ":".join(tokens[0:5])
		if (not merged_snps_hash.get(mykey)):
			#add to hash
			merged_snps_hash[mykey] = strain
		else:
			#print "remove not unique snp: " + mykey
			dupe = dupe + 1
			del merged_snps_hash[mykey]
                	merged_snps_hash[mykey] = strain
	print "collect this statistic and pass to controller removed dupe snps: " + str(dupe)
	file2.close()
	#chr start end sample type notes
	numrecords = len(merged_snps_hash)
	outhash["SNPDupes"] = str(dupe)
	outhash["SNPs"] = str(numrecords)
	processed = 0
	print "processing merging of %i genomic variations" % numrecords	
	for key in merged_snps_hash:
		ktk = key.split(":")		
		mutation_type = getMutationType(merged_snps_hash[key].split(":")[1],ktk[0],ktk[1],ktk[2])	
		if (type_hash.get(mutation_type)):
			mutation_type = type_hash.get(mutation_type)	
		merged_file.write("\t".join(ktk[0:3]) + "\t" + crossing + "_" + merged_snps_hash[key].split(":")[0] + "_SNPMerged\t" + mutation_type + "\t" +  "->".join(ktk[3:5]) + "\n")  
		processed = processed + 1
	merged_file.close() 
	print "Done merge %s" % time.strftime("%c")
	return mutation_repo + "/" + crossing + "_" + pathspl[8] + "_" + pathspl2[8] +  "_merged.mut";

if __name__ == "__main__":
        #Read in config
        #config = ConfigParser.RawConfigParser()
        #config.read('./cgs_lims.config')
	arglen = len(sys.argv)
        #print arglen
	if (arglen < 3):
		print "Proper usage is /tools/bin/python cgs_suite.py strain2bc YSEQ10 outfolder or /tools/bin/python cgs_suite.py cgs_submit /data/mutation/path"
	merge_snps(sys.argv[1], sys.argv[2], sys.argv[3])

